/********************************************/
// Compute transition for MIT shocks..
// this file computes the transitional dynamics after a temporary unexpected shock affecting the self-employment assistance.


#include "PM_method/basculefun.cpp"


#if indexPM == 1
void PM(int T) {    // T is the number of periods, the higher it is, the better is the precision
    
    printf("STARTING Transitional Dynamics...\n");
    
    // INITIALIZATION
    int i, j, numberAGG; // index
    double *pricepath, *nowpricepath, *endprice, *startingprice; // pricepath
    double GOODNESS, critTP;
    double *genmoments, taxout, optimaltaxSEA, capitalout, laborout, fracJJ, fracU, meanwagelike, defaultrate, prodtotal, probleaveU, newENTq; // statistics
    
    
    // FILE
    FILE *pricefile;
    
    // ALLOCATION //
    numberAGG       = 11; // here you specify how any price to follow.
#define inxPrice(indext, indexprice) ((numberAGG)*(indext)+(indexprice))
    
    
    /* PRICEPATH */
    pricepath       = (double *) calloc((T*numberAGG), sizeof(double));
    nowpricepath    = (double *) calloc((T*numberAGG), sizeof(double));
    endprice        = (double *) calloc((numberAGG), sizeof(double));
    startingprice   = (double *) calloc((numberAGG), sizeof(double));
    
    
    startingprice[0] = 0.016510391666088;
    startingprice[1] = 0.253706723228092;
    startingprice[2] = 0.009027144467390;
    endprice[0] = 0.016510391666088;             // Interest rate
    endprice[1] = 0.253706723228092;             // Wage rate
    endprice[2] = 0.009027144467390;             // Tax rate

    // SEA
#if OPT_policy_exp == 1
    endprice[0] = 0.016489214609214;             // Interest rate
    endprice[2] = 0.009314127952320;             // Tax rate
    endprice[1] = 0.253781277852715;             // Wage rate
#endif
    // SEA LS
#if OPT_policy_exp == 2
    endprice[0] = 0.016496051740501;             // Interest rate
    endprice[2] = 0.009294948603129;         // Tax rate
    endprice[1] = 0.253757200218391;             // Wage rate
#endif
    // SEA NB
#if OPT_policy_exp == 3
    endprice[0] = 0.016487125182303;             // Interest rate
    endprice[2] = 0.009425010811369;         // Tax rate
    endprice[1] = 0.253788637357542;             // Wage rate
#endif
    // compensation
#if OPT_policy_exp == 5
    endprice[0] = 0.016489989167641;             // Interest rate
    endprice[2] = 0.009290515499587;         // Tax rate
    endprice[1] = 0.253778549820850;             // Wage rate
#endif
    
    /* MOMENTS */
    genmoments = (double *) calloc((nbmoments), sizeof(double));
    
    // GUESS THE PRICE PATH //
    for(i = 0; i < numberAGG; i++) {
        pricepath[inxPrice(0,i)] = startingprice[i];
        pricepath[inxPrice((T-1),i)] = endprice[i];
    }
    
    for(i = 0; i < numberAGG; i++) {
        for(j = 1; j < (T-1); j++) {
            pricepath[inxPrice(j,i)] = startingprice[i] + j*(endprice[i] - startingprice[i])/(T-1);
            //printf("%f", pricepath[inxPrice(j,i)]); getchar();
        }
    }
    
    /******* SAVE IN FILE - PRICEPATH *******/
    pricefile=fopen(pricepathfile, "w");
    setbuf ( pricefile , NULL );
    fprintf(pricefile,"%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","R", "meanwage", "T", "K", "L", "E", "U", "SEA");
    for(i = 0; i < (T); i++) {
        for(j = 0; j < numberAGG; j++) {
            fprintf(pricefile,"%20.15f\t",pricepath[inxPrice(i, j)]);
        }
        fprintf(pricefile,"\n");
    }
    fclose(pricefile);
    
    
    
#include "PM_method/creatingvalues_CEV.cpp"
#include "PM_method/creatingvalues_PM.cpp"
    
    
    
    
    /******************************************/
    /*** PROJECTED METHOD START REALLY HERE ***/
    /******************************************/
    
    // SAVE THE LAST PERIOD VALUE FUNCTION //
    // E
    bascule2E(valueJET,valueJE,(T-1));
    bascule2E(valueJEiT,valueJEi,(T-1));
    bascule2E(valueJNET,valueJNE,(T-1));
    bascule2E(valueJNEiT,valueJNEi,(T-1));
    // W
    bascule2W(valueWWNET, valueWWNE,(T-1));
    bascule2W(valueWWET, valueWWE,(T-1));
    // U
    bascule2(valueUiNET, valueUiNE,(T-1));
    bascule2(valueUnNET, valueUnNE,(T-1));
    bascule2(valueUiET, valueUiE,(T-1));
    bascule2(valueUnET, valueUnE,(T-1));
    
    
    // SET FIRST PERIOD DISTRIBUTION //
    // E
    bascule2E(distJE_0, distJE, 0);
    bascule2E(distJNE_0, distJNE, 0);
    bascule2E(distJNEi_0, distJNEi, 0);
    bascule2E(distJEi_0, distJEi, 0);
    // W
    bascule2W(distWWE_0, distWWE, 0);
    bascule2W(distWWNE_0, distWWNE, 0);
    // U
    bascule2(distUiNE_0, distUiNE, 0);
    bascule2(distUiE_0, distUiE, 0);
    bascule2(distUnNE_0, distUnNE, 0);
    bascule2(distUnE_0, distUnE, 0);
    
    
    
    /***** STARTING LOOP OVER PRICES *****/
    
    critTP = 1;
    while(critTP > epsilonTP) {
        
        // Plug value in T into the vector corresponding in value //
        bascule(valueJET, temp_valueJE, ifulldimE);
        bascule(valueJNET, temp_valueJNE, ifulldimE);
#if SEA == 1
        bascule(valueJEiT, temp_valueJEi, ifulldimE);
        bascule(valueJNEiT, temp_valueJNEi, ifulldimE);
#endif
        bascule(valueWWNET, temp_valueWWNE, ifulldimW);
        bascule(valueWWET, temp_valueWWE, ifulldimW);
        bascule(valueUiNET, temp_valueUiNE, ifulldim);
        bascule(valueUnNET, temp_valueUnNE, ifulldim);
        bascule(valueUiET, temp_valueUiE, ifulldim);
        bascule(valueUnET, temp_valueUnE, ifulldim);
        
        // SET FIRST PERIOD DISTRIBUTION //
        // E
        bascule(distJE_0, temp_distJE, ifulldimE);
        bascule(distJNE_0, temp_distJNE, ifulldimE);
        bascule(distJNEi_0, temp_distJNEi, ifulldimE);
        bascule(distJEi_0, temp_distJEi, ifulldimE);
        // W
        bascule(distWWE_0, temp_distWWE, ifulldimW);
        bascule(distWWNE_0, temp_distWWNE, ifulldimW);
        // U
        bascule(distUiNE_0, temp_distUiNE, ifulldim);
        bascule(distUiE_0, temp_distUiE, ifulldim);
        bascule(distUnNE_0, temp_distUnNE, ifulldim);
        bascule(distUnE_0, temp_distUnE, ifulldim);
        
        //printf("\n TEST: %20.15f %20.15f ", distWWNE[inxTPW(0, 0, 0, 0)], distWWNE_0[inxW(0, 0, 0)]); getchar();
        
        
        
        // VERIFICATION DISTRIBUTION //
        double distval;
        
        distval = 0.0;
        for(int tgrid = 0; tgrid<maxtheta; tgrid++){
            for(int igrid=0;igrid<maxgrid;igrid++){
                distval+=temp_distUiE[inx(igrid, tgrid)];
                distval+=temp_distUnE[inx(igrid, tgrid)];
                distval+=temp_distUiNE[inx(igrid, tgrid)];
                distval+=temp_distUnNE[inx(igrid, tgrid)];
                
                for(int ygrid=0; ygrid<maxyprod; ygrid++) {
                    distval+=temp_distWWNE[inxW(igrid, tgrid, ygrid)];
                    distval+=temp_distWWE[inxW(igrid, tgrid, ygrid)];
                }
                
                for(int egrid=0; egrid<maxfirmtype; egrid++) {
                    distval+=temp_distJE[inxE(igrid, tgrid, egrid)];
                    distval+=temp_distJNE[inxE(igrid, tgrid, egrid)];
                }
            }
        }
        
        if(distval < 0.9999 || distval > 1.0001) {
            printf("ERROR: %20.15f", distval); getchar();
        }
        
        
        //printf("DIST: %20.15f %20.15f %20.15f %20.15f %20.15f %20.15f %20.15f", temp_distJNE[inxE((maxgrid-1), (maxtheta-1), 0)], temp_distJNE[inxE((maxgrid-1), (maxtheta-1), 1)], temp_distJNE[inxE((maxgrid-1), (maxtheta-1), 2)], temp_distJNE[inxE((maxgrid-1), (maxtheta-1), 3)], temp_distJNE[inxE((maxgrid-1), (maxtheta-1), 4)], temp_distJNE[inxE((maxgrid-1), (maxtheta-1), 5)], temp_distJNE[inxE((maxgrid-1), (maxtheta-1), 6)]);getchar();
        
        
        //////////////////////////////////////////////////////
        // VERIF VALUE FUNCTIONS //
        
        
        /**** COMPUTE POLICY FUNCTIONS AND STORE IT ****/
        
        printf("POLICY VALUE STARTING FOR PM \n");
        
        // i = 0 correspond to i = T-1
        for(int ii = 1; ii < T;  ii++) {
            
            // Set prices right now here //
            rstar       = pricepath[inxPrice((T-1-ii), 0)];
            wstar       = pricepath[inxPrice((T-1-ii), 1)];
            taxrate     = pricepath[inxPrice((T-1-ii), 2)];
            
            printf("PERIOD: %d, R: %f, T: %f, Wage: %f\n", (T-1-ii),  rstar, taxrate, wstar);
            
            /*** VFI FUNCTION ***/
#if NOPOLICY == 1
            VFI(temp_valueUiNE,  temp_saveUiNE, temp_C_valueUiNE, temp_searchUiNE_W, temp_searchUiNE_J,
                temp_valueUiE,   temp_saveUiE,  temp_C_valueUiE,  temp_searchUiE_W, temp_searchUiE_J,
                temp_valueUnNE,  temp_saveUnNE, temp_C_valueUnNE, temp_searchUnNE_W, temp_searchUnNE_J,
                temp_valueUnE,   temp_saveUnE,  temp_C_valueUnE,  temp_searchUnE_W, temp_searchUnE_J,
                temp_valueWWNE,  temp_saveWWNE, temp_C_valueWWNE, temp_searchWWNE_J,
                temp_valueWWE,   temp_saveWWE,  temp_C_valueWWE,  temp_searchWWE_J,
                temp_valueJE,    temp_saveJE,   temp_C_valueJE,   temp_searchJE_W, temp_collateralJE,
                temp_valueJNE,   temp_saveJNE,  temp_C_valueJNE,  temp_searchJNE_W, temp_collateralJNE, temp_sloanJNE,  temp_rateJNE, temp_defaultJNE);
#endif
            
#if SEA == 1
            VFI(temp_valueUiNE,  temp_saveUiNE, temp_C_valueUiNE, temp_searchUiNE_W, temp_searchUiNE_J,
                temp_valueUiE,   temp_saveUiE,  temp_C_valueUiE,  temp_searchUiE_W, temp_searchUiE_J,
                temp_valueUnNE,  temp_saveUnNE, temp_C_valueUnNE, temp_searchUnNE_W, temp_searchUnNE_J,
                temp_valueUnE,   temp_saveUnE,  temp_C_valueUnE,  temp_searchUnE_W, temp_searchUnE_J,
                temp_valueWWNE,  temp_saveWWNE, temp_C_valueWWNE, temp_searchWWNE_J,
                temp_valueWWE,   temp_saveWWE,  temp_C_valueWWE,  temp_searchWWE_J,
                temp_valueJE,    temp_saveJE,   temp_C_valueJE,   temp_searchJE_W, temp_collateralJE,
                temp_valueJNE,   temp_saveJNE,  temp_C_valueJNE,  temp_searchJNE_W, temp_collateralJNE, temp_sloanJNE, temp_rateJNE, temp_defaultJNE,
                temp_valueJEi,   temp_saveJEi,  temp_C_valueJEi,  temp_searchJEi_W, temp_collateralJEi,
                temp_valueJNEi,  temp_saveJNEi, temp_C_valueJNEi, temp_searchJNEi_W, temp_collateralJNEi, temp_sloanJNEi, temp_rateJNEi, temp_defaultJNEi);
#endif
            
            
            /// VALUE FUNCTION ///
            bascule2E(temp_valueJE,valueJE,(T-1-ii));
            bascule2E(temp_valueJEi,valueJEi,(T-1-ii));
            bascule2E(temp_valueJNE,valueJNE,(T-1-ii));
            bascule2E(temp_valueJNEi,valueJNEi,(T-1-ii));
            
            bascule2W(temp_valueWWNE,valueWWNE,(T-1-ii));
            bascule2W(temp_valueWWE,valueWWE,(T-1-ii));
            
            bascule2(temp_valueUiNE,valueUiNE,(T-1-ii));
            bascule2(temp_valueUnNE,valueUnNE,(T-1-ii));
            bascule2(temp_valueUiE,valueUiE,(T-1-ii));
            bascule2(temp_valueUnE,valueUnE,(T-1-ii));
            
            /// POLICY FUNCTION ///
            bascule2E(temp_saveJE,saveJE,(T-1-ii));
            bascule2E(temp_searchJE_W,searchJE_W,(T-1-ii));
            bascule2E(temp_saveJEi,saveJEi,(T-1-ii));
            bascule2E(temp_searchJEi_W,searchJEi_W,(T-1-ii));
            bascule2E(temp_saveJNE,saveJNE,(T-1-ii));
            bascule2E(temp_searchJNE_W,searchJNE_W,(T-1-ii));
            bascule2E(temp_saveJNEi,saveJNEi,(T-1-ii));
            bascule2E(temp_searchJNEi_W,searchJNEi_W,(T-1-ii));
            bascule2E(temp_defaultJNE,defaultJNE,(T-1-ii));
            bascule2E(temp_rateJNE,rateJNE,(T-1-ii));
            bascule2E(temp_defaultJNEi,defaultJNEi,(T-1-ii));
            bascule2E(temp_rateJNEi,rateJNEi,(T-1-ii));
            bascule2E(temp_sloanJNE,sloanJNE,(T-1-ii));
            bascule2E(temp_collateralJNE,collateralJNE,(T-1-ii));
            bascule2E(temp_collateralJE,collateralJE,(T-1-ii));
            bascule2E(temp_sloanJNEi,sloanJNEi,(T-1-ii));
            bascule2E(temp_collateralJNEi,collateralJNEi,(T-1-ii));
            bascule2E(temp_collateralJEi,collateralJEi,(T-1-ii));
            
            bascule2W(temp_saveWWNE,saveWWNE,(T-1-ii));
            bascule2W(temp_searchWWNE_J,searchWWNE_J,(T-1-ii));
            bascule2W(temp_saveWWE,saveWWE,(T-1-ii));
            bascule2W(temp_searchWWE_J,searchWWE_J,(T-1-ii));
            
            bascule2(temp_saveUiNE,saveUiNE,(T-1-ii));
            bascule2(temp_saveUnNE,saveUnNE,(T-1-ii));
            bascule2(temp_saveUiE,saveUiE,(T-1-ii));
            bascule2(temp_saveUnE,saveUnE,(T-1-ii));
            bascule2(temp_searchUiNE_J,searchUiNE_J,(T-1-ii));
            bascule2(temp_searchUnNE_J,searchUnNE_J,(T-1-ii));
            bascule2(temp_searchUnE_J,searchUnE_J,(T-1-ii));
            bascule2(temp_searchUiE_J,searchUiE_J,(T-1-ii));
            bascule2(temp_searchUiNE_W,searchUiNE_W,(T-1-ii));
            bascule2(temp_searchUnNE_W,searchUnNE_W,(T-1-ii));
            bascule2(temp_searchUnE_W,searchUnE_W,(T-1-ii));
            bascule2(temp_searchUiE_W,searchUiE_W,(T-1-ii));
            
        }
        
        
        
        
        /**** COMPUTE DISTRIBUTION AND STORE IT ****/
        // START SIMULATING STARTING FROM PERIOD 0, where prices are given, but policy have changed.//
        
        printf("\n STARTING DISTRIBUTION...\n");
        
        printf("VALUE RANDOM: %f and %f", valueJE[inxTPE((T-1), 14, 1, 3)], valueJNE[inxTPE((T-1), 14, 1, 3)]);
        
        // j = 0 --> moment where policy have changed, but starting distribution is given //
        for(int jj = 0; jj < (T-1); jj++) {
            
            
            // Set prices right now here //
            rstar       = pricepath[inxPrice((jj), 0)];
            wstar       = pricepath[inxPrice((jj), 1)];
            taxrate     = pricepath[inxPrice((jj), 2)];
            
            
            // Set prices right now here //
            bascule3E(valueJE,temp_valueJE,(jj+1));
            bascule3E(valueJEi,temp_valueJEi,(jj+1));
            bascule3E(valueJNE,temp_valueJNE,(jj+1));
            bascule3E(valueJNEi,temp_valueJNEi,(jj+1));
            
            bascule3W(valueWWNE,temp_valueWWNE,(jj+1));
            bascule3W(valueWWE,temp_valueWWE,(jj+1));
            
            bascule3(valueUiNE,temp_valueUiNE,(jj+1));
            bascule3(valueUnNE,temp_valueUnNE,(jj+1));
            bascule3(valueUiE,temp_valueUiE,(jj+1));
            bascule3(valueUnE,temp_valueUnE,(jj+1));
            
            
            // PRINT POLICY INTO TEMP (this is the current behavior) //
            bascule3E(saveJE,temp_saveJE,jj);
            bascule3E(searchJE_W,temp_searchJE_W,jj);
            bascule3E(saveJEi,temp_saveJEi,jj);
            bascule3E(searchJEi_W,temp_searchJEi_W,jj);
            bascule3E(saveJNE,temp_saveJNE,jj);
            bascule3E(searchJNE_W,temp_searchJNE_W,jj);
            bascule3E(saveJNEi,temp_saveJNEi,jj);
            bascule3E(searchJNEi_W,temp_searchJNEi_W,jj);
            bascule3E(defaultJNE,temp_defaultJNE,jj);
            bascule3E(defaultJNEi,temp_defaultJNEi,jj);
            bascule3E(rateJNE,temp_rateJNE,jj);
            bascule3E(rateJNEi,temp_rateJNEi,jj);
            bascule3E(sloanJNE,temp_sloanJNE,jj);
            bascule3E(collateralJNE,temp_collateralJNE,jj);
            bascule3E(collateralJE,temp_collateralJE,jj);
            bascule3E(sloanJNEi,temp_sloanJNEi,jj);
            bascule3E(collateralJNEi,temp_collateralJNEi,jj);
            bascule3E(collateralJEi,temp_collateralJEi,jj);
            
            bascule3W(saveWWNE,temp_saveWWNE,jj);
            bascule3W(searchWWNE_J,temp_searchWWNE_J,jj);
            bascule3W(saveWWE,temp_saveWWE,jj);
            bascule3W(searchWWE_J,temp_searchWWE_J,jj);
            
            bascule3(saveUiNE,temp_saveUiNE,jj);
            bascule3(saveUiE,temp_saveUiE,jj);
            bascule3(saveUnNE,temp_saveUnNE,jj);
            bascule3(saveUnE,temp_saveUnE,jj);
            
            bascule3(searchUiE_W,temp_searchUiE_W,jj);
            bascule3(searchUiE_J,temp_searchUiE_J,jj);
            bascule3(searchUiNE_W,temp_searchUiNE_W,jj);
            bascule3(searchUiNE_J,temp_searchUiNE_J,jj);
            
            bascule3(searchUnNE_J,temp_searchUnNE_J,jj);
            bascule3(searchUnE_J,temp_searchUnE_J,jj);
            bascule3(searchUnNE_W,temp_searchUnNE_W,jj);
            bascule3(searchUnE_W,temp_searchUnE_W,jj);
            
            
            /*** SIMULATION FUNCTION ***/
#if NOPOLICY == 1
            SIMULATION(temp_valueUiNE, temp_C_valueUiNE,  temp_saveUiNE, temp_searchUiNE_W, temp_searchUiNE_J,
                       temp_valueUiE,  temp_C_valueUiE,   temp_saveUiE, temp_searchUiE_W, temp_searchUiE_J,
                       temp_valueUnNE, temp_C_valueUnNE,  temp_saveUnNE, temp_searchUnNE_W, temp_searchUnNE_J,
                       temp_valueUnE,  temp_C_valueUnE,   temp_saveUnE, temp_searchUnE_W, temp_searchUnE_J,
                       temp_valueWWNE, temp_C_valueWWNE,  temp_saveWWNE, temp_searchWWNE_J,
                       temp_valueWWE,  temp_C_valueWWE,   temp_saveWWE, temp_searchWWE_J,
                       temp_valueJE,   temp_C_valueJE,    temp_saveJE, temp_searchJE_W, temp_collateralJE,
                       temp_valueJNE,  temp_C_valueJNE,   temp_saveJNE, temp_searchJNE_W, temp_collateralJNE, temp_sloanJNE, temp_rateJNE, temp_defaultJNE,
                       &capitalout, &laborout, &taxout, genmoments,
                       temp_distJE, temp_distJNE, temp_distWWNE, temp_distWWE, temp_distUnE, temp_distUnNE, temp_distUiNE, temp_distUiE, temp_distJEnext, temp_distJNEnext, temp_distWWNEnext, temp_distWWEnext, temp_distUnEnext, temp_distUnNEnext, temp_distUiNEnext, temp_distUiEnext,
                       &fracJJ, &fracU, &meanwagelike, &optimaltaxSEA, &defaultrate, &prodtotal, &probleaveU, &newENTq);
#endif
            
#if SEA == 1
            SIMULATION(temp_valueUiNE, temp_C_valueUiNE, temp_saveUiNE, temp_searchUiNE_W, temp_searchUiNE_J,
                       temp_valueUiE,  temp_C_valueUiE, temp_saveUiE, temp_searchUiE_W, temp_searchUiE_J,
                       temp_valueUnNE, temp_C_valueUnNE, temp_saveUnNE, temp_searchUnNE_W, temp_searchUnNE_J,
                       temp_valueUnE,  temp_C_valueUnE, temp_saveUnE, temp_searchUnE_W, temp_searchUnE_J,
                       temp_valueWWNE, temp_C_valueWWNE, temp_saveWWNE, temp_searchWWNE_J,
                       temp_valueWWE,  temp_C_valueWWE, temp_saveWWE, temp_searchWWE_J,
                       temp_valueJE,   temp_C_valueJE, temp_saveJE, temp_searchJE_W, temp_collateralJE,
                       temp_valueJNE,  temp_C_valueJNE, temp_saveJNE, temp_searchJNE_W, temp_collateralJNE, temp_sloanJNE, temp_rateJNE, temp_defaultJNE,
                       temp_valueJEi,  temp_C_valueJEi, temp_saveJEi, temp_searchJEi_W, temp_collateralJEi,
                       temp_valueJNEi, temp_C_valueJNEi, temp_saveJNEi, temp_searchJNEi_W, temp_collateralJNEi, temp_sloanJNEi, temp_rateJNEi, temp_defaultJNEi,
                       &capitalout, &laborout, &taxout, &optimaltaxSEA, genmoments, temp_distJE, temp_distJNE, temp_distJEi, temp_distJNEi, temp_distWWNE, temp_distWWE, temp_distUnE, temp_distUnNE, temp_distUiNE, temp_distUiE, temp_distJEnext, temp_distJNEnext, temp_distJEinext, temp_distJNEinext, temp_distWWNEnext, temp_distWWEnext, temp_distUnEnext, temp_distUnNEnext, temp_distUiNEnext, temp_distUiEnext, &fracJJ, &fracU, &meanwagelike, &defaultrate, &prodtotal, &probleaveU, &newENTq);
#endif
            
            
            //printf("K=%f L=%f T=%f FracE=%f fracU=%f, meanwage=%f, default=%f Y=%f Cost=%f",capitalout,laborout,taxout,fracJJ,fracU,meanwagelike,defaultrate,prodtotal,optimaltaxSEA);getchar();
            
            nowpricepath[inxPrice((jj), 0)] = alphapar*TFPpar*pow((laborout/capitalout),(1.0-alphapar)) - deltapar - wedgerate;
            nowpricepath[inxPrice((jj), 1)] = (1.0-alphapar)*TFPpar*pow((capitalout/laborout),(alphapar));
            nowpricepath[inxPrice((jj), 2)] = taxout;
            nowpricepath[inxPrice((jj), 3)] = capitalout;
            nowpricepath[inxPrice((jj), 4)] = 1.0-fracJJ-fracU; //laborout;
            nowpricepath[inxPrice((jj), 5)] = fracJJ;
            nowpricepath[inxPrice((jj), 6)] = fracU;
            nowpricepath[inxPrice((jj), 7)] = 100*(defaultrate);
            nowpricepath[inxPrice((jj), 8)] = prodtotal;
            nowpricepath[inxPrice((jj), 9)] = probleaveU;
            nowpricepath[inxPrice((jj), 10)] = newENTq;
            
            // E
            bascule(temp_distJEnext, temp_distJE, ifulldimE);
            bascule(temp_distJEinext, temp_distJEi, ifulldimE);
            bascule(temp_distJNEnext, temp_distJNE, ifulldimE);
            bascule(temp_distJNEinext, temp_distJNEi, ifulldimE);
            
            bascule(temp_distWWEnext, temp_distWWE, ifulldimW);
            bascule(temp_distWWNEnext, temp_distWWNE, ifulldimW);
            
            bascule(temp_distUnEnext, temp_distUnE, ifulldim);
            bascule(temp_distUnNEnext, temp_distUnNE, ifulldim);
            bascule(temp_distUiNEnext, temp_distUiNE, ifulldim);
            bascule(temp_distUiEnext, temp_distUiE, ifulldim);
            
            bascule2E(temp_distJEnext, distJE, (jj+1));
            bascule2E(temp_distJEinext, distJEi, (jj+1));
            bascule2E(temp_distJNEnext, distJNE, (jj+1));
            bascule2E(temp_distJNEinext, distJNEi, (jj+1));
            
            bascule2W(temp_distWWEnext, distWWE, (jj+1));
            bascule2W(temp_distWWNEnext, distWWNE, (jj+1));
            
            bascule2(temp_distUnEnext, distUnE, (jj+1));
            bascule2(temp_distUnNEnext, distUnNE, (jj+1));
            bascule2(temp_distUiNEnext, distUiNE, (jj+1));
            bascule2(temp_distUiEnext, distUiE, (jj+1));
            
            
            printf("PERIOD: %d \n", jj);
        }
        
        
        //****** UPDATE PRICES up to A RELAXATION ******//
        critTP = 0.0;
        //    for(int ii = 1; ii < (T-1); ii++) {
        //        for(int jj = 0; jj < numberAGG; jj++) {
        //            critTP = max(critTP, abs(pricepath[inxPrice(ii,jj)] - nowpricepath[inxPrice(ii,jj)])/pricepath[inxPrice(ii,jj)]);
        //            pricepath[inxPrice(ii,jj)] = relaxTP*pricepath[inxPrice(ii,jj)]+(1.0 - relaxTP)*nowpricepath[inxPrice(ii,jj)];
        //        }
        //    }
        for(int ii = 1; ii < (T-1); ii++) {
            critTP = max(critTP, abs(pricepath[inxPrice(ii,0)] - nowpricepath[inxPrice(ii,0)])/pricepath[inxPrice(ii,0)]);
            pricepath[inxPrice(ii,0)] = relaxTP*pricepath[inxPrice(ii,0)]+(1.0 - relaxTP)*nowpricepath[inxPrice(ii,0)];
            critTP = max(critTP, abs(pricepath[inxPrice(ii,1)] - nowpricepath[inxPrice(ii,1)])/pricepath[inxPrice(ii,1)]);
            pricepath[inxPrice(ii,1)] = relaxTP*pricepath[inxPrice(ii,1)]+(1.0 - relaxTP)*nowpricepath[inxPrice(ii,1)];
            critTP = max(critTP, abs(pricepath[inxPrice(ii,2)] - nowpricepath[inxPrice(ii,2)])/pricepath[inxPrice(ii,2)]);
            pricepath[inxPrice(ii,2)] = relaxTP*pricepath[inxPrice(ii,2)]+(1.0 - relaxTP)*nowpricepath[inxPrice(ii,2)];
            for(int jj = 3; jj < numberAGG; jj++) {
                pricepath[inxPrice(ii,jj)] = nowpricepath[inxPrice(ii,jj)];
            }
        }
        
        for(j = 0; j < T; j++) {
            for(i = 0; i < numberAGG; i++) {
                printf("%f\t", pricepath[inxPrice(j,i)]);
            }
            printf("\n");
        }
        
        
        
        
        
        //****** COMPUTE THE GOODNESS OF FIT FOR F ******//
        GOODNESS = 0.0;
        for(j = 0; j < numberAGG; j++) {
            GOODNESS = max(GOODNESS, abs(pricepath[inxPrice((T-2), j)] - pricepath[inxPrice((T-1), j)]));
        }
        printf("GOODNESS: %f || CRITERE = %F\n", GOODNESS, critTP);
        
        
        
        
        /******* SAVE IN FILE - PRICEPATH *******/
        pricefile=fopen(pricepathfile, "w");
        setbuf ( pricefile , NULL );
        fprintf(pricefile,"%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","R", "W", "T", "K", "L", "E", "U","B", "Y", "F", "N");
        for(i = 0; i < (T); i++) {
            for(j = 0; j < numberAGG; j++) {
                fprintf(pricefile,"%20.15f\t",pricepath[inxPrice(i, j)]);
            }
            fprintf(pricefile,"\n");
        }
        fclose(pricefile);
        
        
        
        
        
        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
        
        
        
        
        /*******************************************/
        /** COMPUTE THE CEV AND VEIL OF IGNORANCE **/
        /*******************************************/
        
        /* The principle is the following, to compute the CEV, you need following Floden (2011)
         ValB = the new SS life-time utility
         ValA = the earlier SS life-time utility
         
         To compute the CEV in the context of a CRRA utility, you have: wU = (valB/valA)^(...)/(...) */
        
        
        double temp_VALA, temp_VALB, temp_VALA_T, temp_VALC_T, temp_VALB_T, VEILESS, VEILT[T-1];
        
        
        /***************************************************/
        //////////// COMPARISON OF STEADY STATES ////////////
        /***************************************************/
        
        double V_SS_E, V_SS_W, V_SS_U, V_SS_E_nonE, V_SS_E_ui, V_SS_E_nonui_select;
        double totE,totW,totU,totE_nonE,totE_ui,totE_nonui_select;
        
        V_SS_E = 0; V_SS_W = 0; V_SS_U = 0; V_SS_E_nonE = 0; V_SS_E_ui = 0; V_SS_E_nonui_select = 0;
        
        // To compute the veil of ignorance utilitarian welfare change //
        temp_VALA = 0.0;
        temp_VALB = 0.0;
        
        for(int tgrid = 0; tgrid < maxtheta ; tgrid++)  // for every state of the world theta
        {
            for(int igrid=0;igrid<maxgrid;igrid++)  // for any point in the asset distribution
            {
                
                for(int ygrid = 0; ygrid<maxyprod; ygrid++) {   // for any point for the matching shock (only for workers)
                    
                    // WORKER NOT EXCLUDED
                    temp_VALA += distWWNE_0[inxW(igrid,tgrid, ygrid)]*valueWWNE0[inxW(igrid,tgrid, ygrid)];  // compute the weighted value of
                    temp_VALA += distWWE_0[inxW(igrid,tgrid, ygrid)]*valueWWE0[inxW(igrid,tgrid, ygrid)];
                    
                    // WORKER EXCLUDED
                    temp_VALB += distWWNE_T[inxW(igrid,tgrid, ygrid)]*valueWWNET[inxW(igrid,tgrid, ygrid)];
                    temp_VALB += distWWE_T[inxW(igrid,tgrid, ygrid)]*valueWWET[inxW(igrid,tgrid, ygrid)];
                    
                    // compute the conditional welfare change between steady-states //
                    // WORKER NOT EXCLUDED
                    EVWWE_SS[inxW(igrid,tgrid, ygrid)] = pow((valueWWET[inxW(igrid,tgrid, ygrid)]-valueWWE0[inxW(igrid,tgrid, ygrid)])/C_WWE0[inxW(igrid,tgrid, ygrid)]+1,(1/(1-sigmapar))) - 1;
                    EVWWNE_SS[inxW(igrid,tgrid, ygrid)] = pow((valueWWNET[inxW(igrid,tgrid, ygrid)]-valueWWNE0[inxW(igrid,tgrid, ygrid)])/C_WWNE0[inxW(igrid,tgrid, ygrid)]+1,(1/(1-sigmapar))) - 1;
                    
                    
                    V_SS_W += (distWWNE_T[inxW(igrid,tgrid, ygrid)]*valueWWNET[inxW(igrid,tgrid, ygrid)] + distWWE_T[inxW(igrid,tgrid, ygrid)]*valueWWET[inxW(igrid,tgrid, ygrid)]);
                    
                    totW += distWWE_T[inxW(igrid,tgrid, ygrid)] + distWWNE_T[inxW(igrid,tgrid, ygrid)];
                    //        // WORKER NOT EXCLUDED
                    //        EVWWE_T[inxW(igrid,tgrid, ygrid)] = 0.0;    // for later, reinitialise
                    //        EVWWNE_T[inxW(igrid,tgrid, ygrid)] = 0.0;
                }
                
                for(int egrid=0;egrid<maxfirmtype;egrid++) {
                    temp_VALA += distJE_0[inxE(igrid,tgrid, egrid)]*valueJE0[inxE(igrid,tgrid, egrid)];
                    temp_VALB += distJE_T[inxE(igrid,tgrid, egrid)]*valueJET[inxE(igrid,tgrid, egrid)];
                    
                    // compute the conditional welfare change //
                    EVJE_SS[inxE(igrid,tgrid, egrid)] =  pow((valueJET[inxE(igrid,tgrid, egrid)]-valueJE0[inxE(igrid,tgrid, egrid)])/C_JE0[inxE(igrid,tgrid, egrid)]+1,(1/(1-sigmapar))) - 1;
                    
                    temp_VALA += distJNE_0[inxE(igrid,tgrid, egrid)]*valueJNE0[inxE(igrid,tgrid, egrid)];
                    temp_VALB += distJNE_T[inxE(igrid,tgrid, egrid)]*valueJNET[inxE(igrid,tgrid, egrid)];
                    
                    // compute the conditional welfare change at the steady state //
                    EVJNE_SS[inxE(igrid,tgrid, egrid)] =  pow((valueJNET[inxE(igrid,tgrid, egrid)]-valueJNE0[inxE(igrid,tgrid, egrid)])/C_JNE0[inxE(igrid,tgrid, egrid)]+1,(1/(1-sigmapar))) - 1;
                    
                    V_SS_E += distJE_T[inxE(igrid,tgrid, egrid)]*valueJET[inxE(igrid,tgrid, egrid)] + distJNE_T[inxE(igrid,tgrid, egrid)]*valueJNET[inxE(igrid,tgrid, egrid)];
                    
                    V_SS_E_nonE += distJE_T[inxE(igrid,tgrid, egrid)]*valueJET[inxE(igrid,tgrid, egrid)] + distJNE_T[inxE(igrid,tgrid, egrid)]*valueJNET[inxE(igrid,tgrid, egrid)];
                    
                    totE += distJE_T[inxE(igrid,tgrid, egrid)] + distJNE_T[inxE(igrid,tgrid, egrid)];
                    totE_nonE += distJE_T[inxE(igrid,tgrid, egrid)] + distJNE_T[inxE(igrid,tgrid, egrid)];
                    //        EVJE_T[inxE(igrid,tgrid, egrid)] = 0.0;     // for later, reinitialise
                    //        EVJNE_T[inxE(igrid,tgrid, egrid)] = 0.0;
                    
                    
#if SEA == 1
                    temp_VALA += distJEi_0[inxE(igrid,tgrid, egrid)]*valueJE0[inxE(igrid,tgrid, egrid)]; //(does not exist)
                    temp_VALB += distJEi_T[inxE(igrid,tgrid, egrid)]*valueJEiT[inxE(igrid,tgrid, egrid)];
                    
                    // compute the conditional welfare change //
                    EVJEi_SS[inxE(igrid,tgrid, egrid)] =  pow((valueJEiT[inxE(igrid,tgrid, egrid)]-valueJE0[inxE(igrid,tgrid, egrid)])/C_JE0[inxE(igrid,tgrid, egrid)]+1,(1/(1-sigmapar))) - 1;
                    
                    temp_VALA += distJNEi_0[inxE(igrid,tgrid, egrid)]*valueJNE0[inxE(igrid,tgrid, egrid)]; // does not exist yes ok
                    temp_VALB += distJNEi_T[inxE(igrid,tgrid, egrid)]*valueJNEiT[inxE(igrid,tgrid, egrid)];
                    
                    // compute the conditional welfare change at the steady-state //
                    EVJNEi_SS[inxE(igrid,tgrid, egrid)] =  pow((valueJNEiT[inxE(igrid,tgrid, egrid)]-valueJNE0[inxE(igrid,tgrid, egrid)])/C_JNE0[inxE(igrid,tgrid, egrid)]+1,(1/(1-sigmapar))) - 1; // we compare with a situation with no insurance.
                    
                    V_SS_E += distJEi_T[inxE(igrid,tgrid, egrid)]*valueJEiT[inxE(igrid,tgrid, egrid)] + distJNEi_T[inxE(igrid,tgrid, egrid)]*valueJNEiT[inxE(igrid,tgrid, egrid)];
                    
                    V_SS_E_ui += distJEi_T[inxE(igrid,tgrid, egrid)]*valueJEiT[inxE(igrid,tgrid, egrid)] + distJNEi_T[inxE(igrid,tgrid, egrid)]*valueJNEiT[inxE(igrid,tgrid, egrid)];
                    
                    V_SS_E_nonui_select += distJEi_T[inxE(igrid,tgrid, egrid)]*valueJE0[inxE(igrid,tgrid, egrid)] + distJNEi_T[inxE(igrid,tgrid, egrid)]*valueJNE0[inxE(igrid,tgrid, egrid)];
                    
                    totE += distJEi_T[inxE(igrid,tgrid, egrid)] + distJNEi_T[inxE(igrid,tgrid, egrid)];
                    totE_nonui_select += distJEi_T[inxE(igrid,tgrid, egrid)] + distJNEi_T[inxE(igrid,tgrid, egrid)];
                    totE_ui += distJEi_T[inxE(igrid,tgrid, egrid)] + distJNEi_T[inxE(igrid,tgrid, egrid)];
                    
                    //        EVJEi_T[inxE(igrid,tgrid, egrid)] = 0.0;      // for later, reinitialise
                    //        EVJNEi_T[inxE(igrid,tgrid, egrid)] = 0.0;
#endif
                }
                
                // this is completely OK.
                temp_VALA += distUiNE_0[inx(igrid,tgrid)]*valueUiNE0[inx(igrid,tgrid)];
                temp_VALA += distUiE_0[inx(igrid,tgrid)]*valueUiE0[inx(igrid,tgrid)];
                temp_VALA += distUnNE_0[inx(igrid,tgrid)]*valueUnNE0[inx(igrid,tgrid)];
                temp_VALA += distUnE_0[inx(igrid,tgrid)]*valueUnE0[inx(igrid,tgrid)];
                
                // this is completely OK.
                temp_VALB += distUiNE_T[inx(igrid,tgrid)]*valueUiNET[inx(igrid,tgrid)];
                temp_VALB += distUiE_T[inx(igrid,tgrid)]*valueUiET[inx(igrid,tgrid)];
                temp_VALB += distUnNE_T[inx(igrid,tgrid)]*valueUnNET[inx(igrid,tgrid)];
                temp_VALB += distUnE_T[inx(igrid,tgrid)]*valueUnET[inx(igrid,tgrid)];
                
                // compute the conditional welfare change //
                EVUiE_SS[inx(igrid,tgrid)] = pow((valueUiET[inx(igrid,tgrid)]-valueUiE0[inx(igrid,tgrid)])/C_UiE0[inx(igrid,tgrid)]+1,(1/(1-sigmapar))) - 1;
                EVUnE_SS[inx(igrid,tgrid)] = pow((valueUnET[inx(igrid,tgrid)]-valueUnE0[inx(igrid,tgrid)])/C_UnE0[inx(igrid,tgrid)]+1,(1/(1-sigmapar))) - 1;
                EVUiNE_SS[inx(igrid,tgrid)] = pow((valueUiNET[inx(igrid,tgrid)]-valueUiNE0[inx(igrid,tgrid)])/C_UiNE0[inx(igrid,tgrid)]+1,(1/(1-sigmapar))) - 1;
                EVUnNE_SS[inx(igrid,tgrid)] = pow((valueUnNET[inx(igrid,tgrid)]-valueUnNE0[inx(igrid,tgrid)])/C_UnNE0[inx(igrid,tgrid)]+1,(1/(1-sigmapar))) - 1;
                
                V_SS_U += distUiNE_T[inx(igrid,tgrid)]*valueUiNET[inx(igrid,tgrid)] + distUiE_T[inx(igrid,tgrid)]*valueUiET[inx(igrid,tgrid)] + distUnNE_T[inx(igrid,tgrid)]*valueUnNET[inx(igrid,tgrid)] + distUnE_T[inx(igrid,tgrid)]*valueUnET[inx(igrid,tgrid)];
                
                totU += distUiNE_T[inx(igrid,tgrid)] + distUiE_T[inx(igrid,tgrid)] + distUnNE_T[inx(igrid,tgrid)] + distUnE_T[inx(igrid,tgrid)];
                
                //    EVUiE_T[inx(igrid,tgrid)] = 0.0;
                //    EVUnE_T[inx(igrid,tgrid)] = 0.0;
                //    EVUiNE_T[inx(igrid,tgrid)] = 0.0;
                //    EVUnNE_T[inx(igrid,tgrid)] = 0.0;
                
            }
        }
        
        
        
        
        /***************************************************/
        //////// ALONG THE TRANSITIONAL PATH (2nd v) ////////
        /***************************************************/
        
        // objective: compute the welfare change of each guy at the time of the reform, for the three occupation. The reform is not retroactive, so it will be operating only next period (the tax does not rise directly //
        // we need to compute the following: for each individual, the cev will be equal to the (Value under the reform at time 1/Value with no reform at time 1)^(1/(1-\sigma)) - 1
        
        double CEV_E_TOT[maxgrid],CEV_E[maxgrid*maxtheta],CEV_E_TOT_dist[maxgrid],CEV_E_dist[maxgrid*maxtheta];
        double CEV_W_TOT[maxgrid],CEV_W[maxgrid*maxtheta],CEV_W_TOT_dist[maxgrid],CEV_W_dist[maxgrid*maxtheta];
        double CEV_U_TOT[maxgrid],CEV_U[maxgrid*maxtheta],CEV_U_TOT_dist[maxgrid],CEV_U_dist[maxgrid*maxtheta];
        
        for(int igrid=0;igrid<maxgrid;igrid++)
        {
            
            CEV_E_TOT[igrid] = 0;
            CEV_W_TOT[igrid] = 0;
            CEV_U_TOT[igrid] = 0;
            CEV_E_TOT_dist[igrid] = 0;
            CEV_W_TOT_dist[igrid] = 0;
            CEV_U_TOT_dist[igrid] = 0;
            
            for(int tgrid = 0; tgrid < maxtheta ; tgrid++)
            {
                CEV_E[inx(igrid,tgrid)] = 0;
                CEV_W[inx(igrid,tgrid)] = 0;
                CEV_U[inx(igrid,tgrid)] = 0;
                CEV_E_dist[inx(igrid,tgrid)] = 0;
                CEV_W_dist[inx(igrid,tgrid)] = 0;
                CEV_U_dist[inx(igrid,tgrid)] = 0;
            }
        }
        
        
        temp_VALA_T = 0.0;
        temp_VALC_T = 0.0;
        temp_VALB_T = 0.0;
        for(int igrid=0;igrid<maxgrid;igrid++)
        {
            for(int tgrid = 0; tgrid < maxtheta ; tgrid++)
            {
                for(int egrid=0;egrid<maxfirmtype;egrid++)
                {
                    
                    temp_VALA_T += distJE_0[inxE(igrid,tgrid,egrid)]*valueJE0[inxE(igrid,tgrid,egrid)];
                    temp_VALC_T += distJE_0[inxE(igrid,tgrid,egrid)]*C_JE0[inxE(igrid,tgrid,egrid)];
                    temp_VALB_T += distJE_0[inxE(igrid,tgrid,egrid)]*valueJE[inxTPE(0,igrid,tgrid,egrid)];
                    
                    // compute the conditional welfare change //
                    EVJE_T[inxE(igrid,tgrid,egrid)] = pow((valueJE[inxTPE(0,igrid,tgrid,egrid)]-valueJE0[inxE(igrid,tgrid,egrid)])/C_JE0[inxE(igrid,tgrid,egrid)]+1,(1/(1-sigmapar))) - 1;
                    
                    
                    temp_VALA_T += distJNE_0[inxE(igrid,tgrid,egrid)]*valueJNE0[inxE(igrid,tgrid,egrid)];
                    temp_VALC_T += distJNE_0[inxE(igrid,tgrid,egrid)]*C_JNE0[inxE(igrid,tgrid,egrid)];
                    temp_VALB_T += distJNE_0[inxE(igrid,tgrid,egrid)]*valueJNE[inxTPE(0,igrid,tgrid,egrid)];
                    
                    // compute the conditional welfare change //
                    EVJNE_T[inxE(igrid,tgrid, egrid)] = pow((valueJNE[inxTPE(0,igrid,tgrid,egrid)]-valueJNE0[inxE(igrid,tgrid,egrid)])/C_JNE0[inxE(igrid,tgrid,egrid)]+1,(1/(1-sigmapar))) - 1;
                    
                    CEV_E_TOT[igrid] += (EVJNE_T[inxE(igrid,tgrid, egrid)]*distJNE_0[inxE(igrid,tgrid,egrid)] + EVJE_T[inxE(igrid,tgrid, egrid)]*distJE_0[inxE(igrid,tgrid,egrid)])/(distJNE_0[inxE(igrid,tgrid,egrid)]+distJE_0[inxE(igrid,tgrid,egrid)]);
                    CEV_E[inx(igrid,tgrid)] += (EVJNE_T[inxE(igrid,tgrid,egrid)]*distJNE_0[inxE(igrid,tgrid,egrid)] + EVJE_T[inxE(igrid,tgrid, egrid)]*distJE_0[inxE(igrid,tgrid,egrid)])/(distJNE_0[inxE(igrid,tgrid,egrid)]+distJE_0[inxE(igrid,tgrid,egrid)]);
                    
                    CEV_E_TOT_dist[igrid] += (distJNE_0[inxE(igrid,tgrid,egrid)]+distJE_0[inxE(igrid,tgrid,egrid)]);
                    CEV_E_dist[inx(igrid,tgrid)] += (distJNE_0[inxE(igrid,tgrid,egrid)]+distJE_0[inxE(igrid,tgrid,egrid)]);
                    //            #if SEA == 1
                    //            temp_VALA_T += distJEi[inxTPE((ii-1), igrid, tgrid, egrid)]*valueJEi[inxTPE((ii-1), igrid, tgrid, egrid)];
                    //            temp_VALB_T += distJEi[inxTPE((ii), igrid, tgrid, egrid)]*valueJEi[inxTPE((ii), igrid, tgrid, egrid)];
                    //
                    //            // compute the conditional welfare change //
                    //            EVJEi_T[inxE(igrid,tgrid, egrid)] += pow((valueJEi[inxTPE((ii), igrid, tgrid, egrid)]/valueJEi[inxTPE((ii-1), igrid, tgrid, egrid)]),(1/(1-sigmapar))) - 1;
                    //            temp_VALA_T += distJNEi[inxTPE((ii-1), igrid, tgrid, egrid)]*valueJNEi[inxTPE((ii-1), igrid, tgrid, egrid)];
                    //            temp_VALB_T += distJNEi[inxTPE((ii), igrid, tgrid, egrid)]*valueJNEi[inxTPE((ii), igrid, tgrid, egrid)];
                    //
                    //            // compute the conditional welfare change //
                    //            EVJNEi_T[inxE(igrid,tgrid, egrid)] += pow((valueJNEi[inxTPE((ii), igrid, tgrid, egrid)]/valueJNEi[inxTPE((ii-1), igrid, tgrid, egrid)]),(1/(1-sigmapar))) - 1;
                    //            #endif
                }
                
                for(int ygrid = 0; ygrid <maxyprod; ygrid++) {
                    
                    temp_VALA_T += distWWNE_0[inxW(igrid,tgrid,ygrid)]*valueWWNE0[inxW(igrid,tgrid,ygrid)];
                    temp_VALA_T += distWWE_0[inxW(igrid,tgrid,ygrid)]*valueWWE0[inxW(igrid,tgrid,ygrid)];
                    temp_VALB_T += distWWNE_0[inxW(igrid,tgrid,ygrid)]*valueWWNE[inxTPW(0,igrid,tgrid,ygrid)];
                    temp_VALB_T += distWWE_0[inxW(igrid,tgrid,ygrid)]*valueWWE[inxTPW(0,igrid,tgrid,ygrid)];

                    temp_VALC_T += distWWNE_0[inxW(igrid,tgrid,ygrid)]*C_WWNE0[inxW(igrid,tgrid,ygrid)];
                    temp_VALC_T += distWWE_0[inxW(igrid,tgrid,ygrid)]*C_WWE0[inxW(igrid,tgrid,ygrid)];

                    // compute the conditional welfare change //
                    EVWWE_T[inxW(igrid,tgrid, ygrid)] = pow((valueWWE[inxTPW(0,igrid,tgrid,ygrid)]-valueWWE0[inxW(igrid,tgrid,ygrid)])/C_WWE0[inxW(igrid,tgrid,ygrid)]+1,(1/(1-sigmapar))) - 1;
                    EVWWNE_T[inxW(igrid,tgrid, ygrid)] = pow((valueWWNE[inxTPW(0,igrid,tgrid,ygrid)]-valueWWNE0[inxW(igrid,tgrid,ygrid)])/C_WWNE0[inxW(igrid,tgrid,ygrid)]+1,(1/(1-sigmapar))) - 1;
                    
                    CEV_W_TOT[igrid] += (EVWWE_T[inxW(igrid,tgrid, ygrid)]*distWWE_0[inxW(igrid,tgrid,ygrid)] + EVWWNE_T[inxW(igrid,tgrid, ygrid)]*distWWNE_0[inxW(igrid,tgrid,ygrid)])/(distWWE_0[inxW(igrid,tgrid,ygrid)] + distWWNE_0[inxW(igrid,tgrid,ygrid)]);
                    CEV_W[inx(igrid,tgrid)] += (EVWWE_T[inxW(igrid,tgrid, ygrid)]*distWWE_0[inxW(igrid,tgrid,ygrid)] + EVWWNE_T[inxW(igrid,tgrid, ygrid)]*distWWNE_0[inxW(igrid,tgrid,ygrid)])/(distWWE_0[inxW(igrid,tgrid,ygrid)] + distWWNE_0[inxW(igrid,tgrid,ygrid)]);
                    
                    CEV_W_TOT_dist[igrid] += distWWE_0[inxW(igrid,tgrid,ygrid)] + distWWNE_0[inxW(igrid,tgrid,ygrid)];
                    CEV_W_dist[inx(igrid,tgrid)] += distWWE_0[inxW(igrid,tgrid,ygrid)] + distWWNE_0[inxW(igrid,tgrid,ygrid)];
                    
                }
                
                
                temp_VALA_T += distUiNE_0[inx(igrid,tgrid)]*valueUiNE0[inx(igrid,tgrid)];
                temp_VALA_T += distUiE_0[inx(igrid,tgrid)]*valueUiE0[inx(igrid,tgrid)];
                temp_VALA_T += distUnNE_0[inx(igrid,tgrid)]*valueUnNE0[inx(igrid,tgrid)];
                temp_VALA_T += distUnE_0[inx(igrid,tgrid)]*valueUnE0[inx(igrid,tgrid)];
                
                temp_VALB_T += distUiNE_0[inx(igrid,tgrid)]*valueUiNE[inxTP(0,igrid,tgrid)];
                temp_VALB_T += distUiE_0[inx(igrid,tgrid)]*valueUiE[inxTP(0,igrid,tgrid)];
                temp_VALB_T += distUnNE_0[inx(igrid,tgrid)]*valueUnNE[inxTP(0,igrid,tgrid)];
                temp_VALB_T += distUnE_0[inx(igrid,tgrid)]*valueUnE[inxTP(0,igrid,tgrid)];

                temp_VALC_T += distUiNE_0[inx(igrid,tgrid)]*C_UiNE0[inx(igrid,tgrid)];
                temp_VALC_T += distUiE_0[inx(igrid,tgrid)]*C_UiE0[inx(igrid,tgrid)];
                temp_VALC_T += distUnNE_0[inx(igrid,tgrid)]*C_UnNE0[inx(igrid,tgrid)];
                temp_VALC_T += distUnE_0[inx(igrid,tgrid)]*C_UnE0[inx(igrid,tgrid)];
                
                // compute the conditional welfare change //
                EVUiE_T[inx(igrid,tgrid)] = pow((valueUiE[inxTP(0,igrid,tgrid)]-valueUiE0[inx(igrid,tgrid)])/C_UiE0[inx(igrid,tgrid)]+1,(1/(1-sigmapar))) - 1;
                EVUnE_T[inx(igrid,tgrid)] = pow((valueUnE[inxTP(0,igrid,tgrid)]-valueUnE0[inx(igrid,tgrid)])/C_UnE0[inx(igrid,tgrid)]+1,(1/(1-sigmapar))) - 1;
                EVUiNE_T[inx(igrid,tgrid)] = pow((valueUiNE[inxTP(0,igrid,tgrid)]-valueUiNE0[inx(igrid,tgrid)])/C_UiNE0[inx(igrid,tgrid)]+1,(1/(1-sigmapar))) - 1;
                EVUnNE_T[inx(igrid,tgrid)] = pow((valueUnNE[inxTP(0,igrid,tgrid)]-valueUnNE0[inx(igrid,tgrid)])/C_UnNE0[inx(igrid,tgrid)]+1,(1/(1-sigmapar))) - 1;
                
                CEV_U_TOT[igrid] += (EVUiE_T[inx(igrid,tgrid)]*distUiE_0[inx(igrid,tgrid)] + EVUiNE_T[inx(igrid,tgrid)]*distUiNE_0[inx(igrid,tgrid)] + EVUnE_T[inx(igrid,tgrid)]*distUnE_0[inx(igrid,tgrid)] + EVUnNE_T[inx(igrid,tgrid)]*distUnNE_0[inx(igrid,tgrid)])/(distUiE_0[inx(igrid,tgrid)] + distUiNE_0[inx(igrid,tgrid)] + distUnE_0[inx(igrid,tgrid)] + distUnNE_0[inx(igrid,tgrid)]);
                CEV_U[inx(igrid,tgrid)] += (EVUiE_T[inx(igrid,tgrid)]*distUiE_0[inx(igrid,tgrid)] + EVUiNE_T[inx(igrid,tgrid)]*distUiNE_0[inx(igrid,tgrid)] + EVUnE_T[inx(igrid,tgrid)]*distUnE_0[inx(igrid,tgrid)] + EVUnNE_T[inx(igrid,tgrid)]*distUnNE_0[inx(igrid,tgrid)])/(distUiE_0[inx(igrid,tgrid)] + distUiNE_0[inx(igrid,tgrid)] + distUnE_0[inx(igrid,tgrid)] + distUnNE_0[inx(igrid,tgrid)]);
                CEV_U_TOT_dist[igrid] += distUiE_0[inx(igrid,tgrid)] + distUiNE_0[inx(igrid,tgrid)] + distUnE_0[inx(igrid,tgrid)] + distUnNE_0[inx(igrid,tgrid)];
                CEV_U_dist[inx(igrid,tgrid)] += distUiE_0[inx(igrid,tgrid)] + distUiNE_0[inx(igrid,tgrid)] + distUnE_0[inx(igrid,tgrid)] + distUnNE_0[inx(igrid,tgrid)];
                
            }
        }
        
        
        // UTILITARIAN COMPARISON //
        VEILESS = pow((temp_VALB/temp_VALA),(1/(1-sigmapar))) - 1;
        VEILT[0] = pow((temp_VALB_T-temp_VALA_T)/temp_VALC_T+1,(1/(1-sigmapar))) - 1; // this should be ok.
        
        
        
        printf("\n VEILESS = %f [valA=%f, valB=%F],  VEILT = %5.10f [valA=%f, valB=%f, valC=%f] \n", VEILESS*100,  temp_VALA, temp_VALB, VEILT[0]*100, temp_VALA_T, temp_VALB_T, temp_VALC_T);
        
        printf("\n TEST U | SS: %f %f %f | TRANSI: %f %f %f \n", EVUiNE_SS[inx(30,1)]*100, valueUiNET[inx(30,1)], valueUiNE0[inx(30,1)], EVUiNE_T[inx(30,1)]*100, valueUiNE[inxTP(0,30,1)], valueUiNE0[inx(30,1)]);
        printf("\n TEST W | SS: %f %f %f | TRANSI: %f %f %f \n", EVWWNE_SS[inxW(30,1,3)]*100, valueWWNET[inxW(30,1,3)], valueWWNE0[inxW(30,1,3)], EVWWNE_T[inxW(30,1,3)]*100, valueWWNE[inxTPW(0,30,1,3)], valueWWNE0[inxW(30,1,3)]);
        
        
        
        
        /***************************************************/
        //////////// ALONG THE TRANSITIONAL PATH ////////////
        /***************************************************/
        
        //for(int ii = 1; ii < (T-1);  ii++)
        //{
        //
        //// REINITIALIZE CRITERION //
        //temp_VALA_T = 0.0;
        //temp_VALB_T = 0.0;
        //
        //    for(int tgrid = 0; tgrid < maxtheta ; tgrid++)
        //    {
        //        for(int igrid=0;igrid<maxgrid;igrid++)
        //        {
        //            for(int egrid=0;egrid<maxfirmtype;egrid++)
        //            {
        //
        //                temp_VALA_T += distJE[inxTPE((ii-1), igrid, tgrid, egrid)]*valueJE[inxTPE((ii-1), igrid, tgrid, egrid)];
        //                temp_VALB_T += distJE[inxTPE((ii), igrid, tgrid, egrid)]*valueJE[inxTPE((ii), igrid, tgrid, egrid)];
        //
        //                // compute the conditional welfare change //
        //                EVJE_T[inxE(igrid,tgrid, egrid)] += pow((valueJE[inxTPE((ii), igrid, tgrid, egrid)]/valueJE[inxTPE((ii-1), igrid, tgrid, egrid)]),(1/(1-sigmapar))) - 1;
        //
        //                temp_VALA_T += distJNE[inxTPE((ii-1), igrid, tgrid, egrid)]*valueJNE[inxTPE((ii-1), igrid, tgrid, egrid)];
        //                temp_VALB_T += distJNE[inxTPE((ii), igrid, tgrid, egrid)]*valueJNE[inxTPE((ii), igrid, tgrid, egrid)];
        //
        //                // compute the conditional welfare change //
        //                EVJNE_T[inxE(igrid,tgrid, egrid)] += pow((valueJNE[inxTPE((ii), igrid, tgrid, egrid)]/valueJNE[inxTPE((ii-1), igrid, tgrid, egrid)]),(1/(1-sigmapar))) - 1;
        //
        //                #if SEA == 1
        //                temp_VALA_T += distJEi[inxTPE((ii-1), igrid, tgrid, egrid)]*valueJEi[inxTPE((ii-1), igrid, tgrid, egrid)];
        //                temp_VALB_T += distJEi[inxTPE((ii), igrid, tgrid, egrid)]*valueJEi[inxTPE((ii), igrid, tgrid, egrid)];
        //
        //                // compute the conditional welfare change //
        //                EVJEi_T[inxE(igrid,tgrid, egrid)] += pow((valueJEi[inxTPE((ii), igrid, tgrid, egrid)]/valueJEi[inxTPE((ii-1), igrid, tgrid, egrid)]),(1/(1-sigmapar))) - 1;
        //                temp_VALA_T += distJNEi[inxTPE((ii-1), igrid, tgrid, egrid)]*valueJNEi[inxTPE((ii-1), igrid, tgrid, egrid)];
        //                temp_VALB_T += distJNEi[inxTPE((ii), igrid, tgrid, egrid)]*valueJNEi[inxTPE((ii), igrid, tgrid, egrid)];
        //
        //                // compute the conditional welfare change //
        //                EVJNEi_T[inxE(igrid,tgrid, egrid)] += pow((valueJNEi[inxTPE((ii), igrid, tgrid, egrid)]/valueJNEi[inxTPE((ii-1), igrid, tgrid, egrid)]),(1/(1-sigmapar))) - 1;
        //                #endif
        //            }
        //
        //            for(int ygrid = 0; ygrid <maxyprod; ygrid++) {
        //                temp_VALA_T += distWWNE[inxTPW((ii-1), igrid, tgrid, ygrid)]*valueWWNE[inxTPW((ii-1), igrid, tgrid, ygrid)];
        //                temp_VALA_T += distWWE[inxTPW((ii-1), igrid, tgrid, ygrid)]*valueWWE[inxTPW((ii-1), igrid, tgrid, ygrid)];
        //                temp_VALB_T += distWWNE[inxTPW((ii), igrid, tgrid, ygrid)]*valueWWNE[inxTPW((ii), igrid, tgrid, ygrid)];
        //                temp_VALB_T += distWWE[inxTPW((ii), igrid, tgrid, ygrid)]*valueWWE[inxTPW((ii), igrid, tgrid, ygrid)];
        //
        //                // compute the conditional welfare change //
        //                EVWWE_T[inxW(igrid,tgrid, ygrid)] += pow((valueWWE[inxTPW((ii), igrid, tgrid, ygrid)]/valueWWE[inxTPW((ii-1), igrid, tgrid, ygrid)]),(1/(1-sigmapar))) - 1;
        //                EVWWNE_T[inxW(igrid,tgrid, ygrid)] += pow((valueWWNE[inxTPW((ii), igrid, tgrid, ygrid)]/valueWWNE[inxTPW((ii-1), igrid, tgrid, ygrid)]),(1/(1-sigmapar))) - 1;
        //            }
        //
        //
        //        temp_VALA_T += distUiNE[inxTP((ii-1), igrid, tgrid)]*valueUiNE[inxTP((ii-1), igrid, tgrid)];
        //        temp_VALA_T += distUiE[inxTP((ii-1), igrid, tgrid)]*valueUiE[inxTP((ii-1), igrid, tgrid)];
        //        temp_VALA_T += distUnNE[inxTP((ii-1), igrid, tgrid)]*valueUnNE[inxTP((ii-1), igrid, tgrid)];
        //        temp_VALA_T += distUnE[inxTP((ii-1), igrid, tgrid)]*valueUnE[inxTP((ii-1), igrid, tgrid)];
        //
        //        temp_VALB_T += distUiNE[inxTP((ii), igrid, tgrid)]*valueUiNE[inxTP((ii), igrid, tgrid)];
        //        temp_VALB_T += distUiE[inxTP((ii), igrid, tgrid)]*valueUiE[inxTP((ii), igrid, tgrid)];
        //        temp_VALB_T += distUnNE[inxTP((ii), igrid, tgrid)]*valueUnNE[inxTP((ii), igrid, tgrid)];
        //        temp_VALB_T += distUnE[inxTP((ii), igrid, tgrid)]*valueUnE[inxTP((ii), igrid, tgrid)];
        //
        //        // compute the conditional welfare change //
        //        EVUiE_T[inx(igrid,tgrid)] += pow((valueUiE[inxTP((ii), igrid, tgrid)]/valueUiE[inxTP((ii-1), igrid, tgrid)]),(1/(1-sigmapar))) - 1;
        //        EVUnE_T[inx(igrid,tgrid)] += pow((valueUnE[inxTP((ii), igrid, tgrid)]/valueUnE[inxTP((ii-1), igrid, tgrid)]),(1/(1-sigmapar))) - 1;
        //        EVUiNE_T[inx(igrid,tgrid)] += pow((valueUiNE[inxTP((ii), igrid, tgrid)]/valueUiNE[inxTP((ii-1), igrid, tgrid)]),(1/(1-sigmapar))) - 1;
        //        EVUnNE_T[inx(igrid,tgrid)] += pow((valueUnNE[inxTP((ii), igrid, tgrid)]/valueUnNE[inxTP((ii-1), igrid, tgrid)]),(1/(1-sigmapar))) - 1;
        //
        //        }
        //    }
        //
        //if((temp_VALB_T/temp_VALA_T) < 0 | temp_VALA_T == 0.0){printf("ERROR CEV: %f", (temp_VALB_T/temp_VALA_T));getchar();}
        //VEILT[ii] = pow((temp_VALB_T/temp_VALA_T),(1/(1-sigmapar))) - 1;
        ////printf("\n VEIL: %20.15f %20.15f %20.15f %20.15f %20.15f %20.15f %20.15f", VEILT[ii], temp_VALB_T, temp_VALA_T, distWWNE[inxTPW((ii-1), 0, 0, 0)], distWWNE[inxTPW((ii), 0, 0, 0)], valueWWNE[inxTPW((ii-1), 0, 0, 0)], valueWWNE[inxTPW((ii), 0, 0, 0)]); getchar();
        //}
        
        /****** END COMPUTATION OF CEV ******/
        
        
        
        
        /************************/
        /** PRINT RESULTS HERE **/
        /************************/
        
        FILE *welfarefile;

        
        // PRINT RESULTS SS VS SS //
        welfarefile=fopen(VEILTfile, "w");
        setbuf ( welfarefile , NULL );
        //fprintf(welfarefile,"%20s\n","VEILT");
        //for(int tt = 0; tt < (T-1) ; tt++)
        //{
        //    fprintf(welfarefile,"%20.15f\n",);
        //}
        fprintf(welfarefile,"\n VEILESS = %f [valA=%f, valB=%F],  VEILT = %5.10f [valA=%f, valB=%f, valC=%f] \n", VEILESS,  temp_VALA, temp_VALB, VEILT[0], temp_VALA_T, temp_VALB_T, temp_VALC_T);
        fprintf(welfarefile,"V_SS_E = %20.15f, V_SS_W = %20.15f, V_SS_U = %20.15f, V_SS_E_nonE = %20.15f, V_SS_E_ui = %20.15f, V_SS_E_nonui_select = %20.15f,\n",V_SS_E, V_SS_W, V_SS_U, V_SS_E_nonE, V_SS_E_ui, V_SS_E_nonui_select);
        fprintf(welfarefile,"totE = %20.15f, totW = %20.15f, totU = %20.15f, totE_nonE = %20.15f, totE_ui = %20.15f, totE_nonui_select = %20.15f,\n", totE,totW,totU,totE_nonE,totE_ui,totE_nonui_select);
        fprintf(welfarefile,"totE = %20.15f, totW = %20.15f, totU = %20.15f, totE_nonE = %20.15f, totE_ui = %20.15f, totE_nonui_select = %20.15f,\n", V_SS_E/totE,V_SS_W/totW,V_SS_U/totU,V_SS_E_nonE/totE_nonE,V_SS_E_ui/totE_ui,V_SS_E_nonui_select/totE_nonui_select);
        fclose(welfarefile);
        
        
        
        
        
        // PRINT RESULTS CEV FOR IGRID //
        welfarefile=fopen(CEV_TOT, "w");
        setbuf ( welfarefile , NULL );
        fprintf(welfarefile,"%20s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","ASSET", "THETA", "distU", "Welf_U", "distW","Welf_W", "distE", "Welf_E");
        for(int tgrid = 0; tgrid < maxtheta ; tgrid++)
        {
            for(int igrid=0;igrid<maxgrid;igrid++)
            {
                fprintf(welfarefile,"%20.15f\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",grid[igrid], tgrid, CEV_U_dist[inx(igrid, tgrid)],CEV_U[inx(igrid,tgrid)], CEV_W_dist[inx(igrid, tgrid)],CEV_W[inx(igrid,tgrid)],CEV_E_dist[inx(igrid, tgrid)],CEV_E[inx(igrid,tgrid)]);
                
            }
        }
        fclose(welfarefile);
        
        
        // PRINT RESULTS CEV FOR IGRID / TGRID //
        welfarefile=fopen(CEV, "w");
        setbuf ( welfarefile , NULL );
        fprintf(welfarefile,"%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","ASSET", "distU", "Welf_U", "distW","Welf_W", "distE", "Welf_E");
        for(int igrid=0;igrid<maxgrid;igrid++)
        {
            fprintf(welfarefile,"%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",grid[igrid], CEV_U_TOT_dist[igrid],CEV_U_TOT[igrid], CEV_W_TOT_dist[igrid],CEV_W_TOT[igrid],CEV_E_TOT_dist[igrid],CEV_E_TOT[igrid]);
            
        }
        fclose(welfarefile);
        
        
        
        
        
        
        // PRINT RESULTS CEV TAKING INTO ACCOUNT TRANSITION FOR ALL //
        // UNEMPLOYED //
        welfarefile=fopen(CEV_U_file, "w");
        setbuf ( welfarefile , NULL );
        fprintf(welfarefile,"%20s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","ASSET", "PROD", "distUiE", "SSUiE", "TUiE","distUiNE", "SSUiNE", "TUiNE","distUnE", "SSUnE", "TUnE","distUnNE", "SSUnNE", "TUnNE");
        for(int tgrid = 0; tgrid < maxtheta ; tgrid++)
        {
            for(int igrid=0;igrid<maxgrid;igrid++)
            {
                fprintf(welfarefile,"%20.15f\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",grid[igrid], tgrid, distUiE_0[inx(igrid, tgrid)],EVUiE_SS[inx(igrid,tgrid)], EVUiE_T[inx(igrid,tgrid)],distUiNE_0[inx(igrid, tgrid)],EVUiNE_SS[inx(igrid,tgrid)], EVUiNE_T[inx(igrid,tgrid)],distUnE_0[inx(igrid, tgrid)],EVUnE_SS[inx(igrid,tgrid)], EVUnE_T[inx(igrid,tgrid)],distUnNE_0[inx(igrid, tgrid)],EVUnNE_SS[inx(igrid,tgrid)], EVUnNE_T[inx(igrid,tgrid)] );
                
            }
        }
        fclose(welfarefile);
        
        // WORKER //
        welfarefile=fopen(CEV_W_file, "w");
        setbuf ( welfarefile , NULL );
        fprintf(welfarefile,"%20s\t%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","ASSET", "PROD", "YPROD", "distWWE", "SSWWE", "TWWE","distWWNE", "SSWWNE", "TWWNE");
        for(int tgrid = 0; tgrid < maxtheta ; tgrid++)
        {
            for(int igrid=0;igrid<maxgrid;igrid++)
            {
                for(int ygrid = 0; ygrid<maxyprod; ygrid++)
                {
                    fprintf(welfarefile,"%20.15f\t%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",grid[igrid], tgrid, ygrid, distWWE_0[inxW(igrid, tgrid, ygrid)],EVWWE_SS[inxW(igrid, tgrid, ygrid)], EVWWE_T[inxW(igrid, tgrid, ygrid)],distWWNE_0[inxW(igrid, tgrid, ygrid)],EVWWNE_SS[inxW(igrid, tgrid, ygrid)], EVWWNE_T[inxW(igrid, tgrid, ygrid)]);
                }
            }
        }
        fclose(welfarefile);
        
        // ENTREPRENEUR //
#if SEA == 1
        welfarefile=fopen(CEV_J, "w");
        setbuf ( welfarefile , NULL );
        fprintf(welfarefile,"%20s\t%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","ASSET","TYPE", "PROD", "distJE", "SSJE", "TJE", "distJEi", "SSJEi", "TJEi", "distJNE", "SSJNE", "TJNE", "distJNEi", "SSJNEi", "TJNEi");
        for(int tgrid = 0; tgrid < maxtheta ; tgrid++)
        {
            for(int igrid=0;igrid<maxgrid;igrid++)
            {
                for(int egrid=0;egrid<maxfirmtype;egrid++)
                {
                    fprintf(welfarefile,"%20.15f\t%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",grid[igrid], egrid, tgrid, distJE_0[inxE(igrid, tgrid,egrid)],EVJE_SS[inxE(igrid, tgrid,egrid)], EVJE_T[inxE(igrid, tgrid,egrid)], distJEi_0[inxE(igrid, tgrid,egrid)],EVJEi_SS[inxE(igrid, tgrid,egrid)], EVJEi_T[inxE(igrid, tgrid,egrid)], distJNE_0[inxE(igrid, tgrid,egrid)],EVJE_SS[inxE(igrid, tgrid,egrid)], EVJNE_T[inxE(igrid, tgrid,egrid)], distJNEi_0[inxE(igrid, tgrid,egrid)],EVJNEi_SS[inxE(igrid, tgrid,egrid)], EVJNEi_T[inxE(igrid, tgrid,egrid)]);
                }
                
            }
        }
        fclose(welfarefile);
#endif
#if SEA == 0
        welfarefile=fopen(CEV_J, "w");
        setbuf ( welfarefile , NULL );
        fprintf(welfarefile,"%20s\t%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","ASSET","TYPE", "PROD", "distJE", "SSJE", "TJE", "distJNE", "SSJNE", "TJNE");
        for(int tgrid = 0; tgrid < maxtheta ; tgrid++)
        {
            for(int igrid=0;igrid<maxgrid;igrid++)
            {
                for(int egrid=0;egrid<maxfirmtype;egrid++)
                {
                    fprintf(welfarefile,"%20.15f\t%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",grid[igrid], egrid, tgrid, distJE_0[inxE(igrid, tgrid,egrid)],EVJE_SS[inxE(igrid, tgrid,egrid)], EVJE_T[inxE(igrid, tgrid,egrid)], distJNE_0[inxE(igrid, tgrid,egrid)],EVJE_SS[inxE(igrid, tgrid,egrid)], EVJNE_T[inxE(igrid, tgrid,egrid)]);
                }
                
            }
        }
        fclose(welfarefile);
#endif
        
        
        
        printf("Crit cvng: %20.15f, %20.15f \n", critTP, VEILESS);
        
    } // end while over TP
    
    
#include "PM_method/freepointers.cpp"
    
}
#endif
